home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dorgbr.f < prev    next >
Text File  |  1996-07-19  |  7KB  |  224 lines
  1.       SUBROUTINE DORGBR( VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO )
  2. *
  3. *  -- LAPACK routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     September 30, 1994
  7. *
  8. *     .. Scalar Arguments ..
  9.       CHARACTER          VECT
  10.       INTEGER            INFO, K, LDA, LWORK, M, N
  11. *     ..
  12. *     .. Array Arguments ..
  13.       DOUBLE PRECISION   A( LDA, * ), TAU( * ), WORK( LWORK )
  14. *     ..
  15. *
  16. *  Purpose
  17. *  =======
  18. *
  19. *  DORGBR generates one of the real orthogonal matrices Q or P**T
  20. *  determined by DGEBRD when reducing a real matrix A to bidiagonal
  21. *  form: A = Q * B * P**T.  Q and P**T are defined as products of
  22. *  elementary reflectors H(i) or G(i) respectively.
  23. *
  24. *  If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q
  25. *  is of order M:
  26. *  if m >= k, Q = H(1) H(2) . . . H(k) and DORGBR returns the first n
  27. *  columns of Q, where m >= n >= k;
  28. *  if m < k, Q = H(1) H(2) . . . H(m-1) and DORGBR returns Q as an
  29. *  M-by-M matrix.
  30. *
  31. *  If VECT = 'P', A is assumed to have been a K-by-N matrix, and P**T
  32. *  is of order N:
  33. *  if k < n, P**T = G(k) . . . G(2) G(1) and DORGBR returns the first m
  34. *  rows of P**T, where n >= m >= k;
  35. *  if k >= n, P**T = G(n-1) . . . G(2) G(1) and DORGBR returns P**T as
  36. *  an N-by-N matrix.
  37. *
  38. *  Arguments
  39. *  =========
  40. *
  41. *  VECT    (input) CHARACTER*1
  42. *          Specifies whether the matrix Q or the matrix P**T is
  43. *          required, as defined in the transformation applied by DGEBRD:
  44. *          = 'Q':  generate Q;
  45. *          = 'P':  generate P**T.
  46. *
  47. *  M       (input) INTEGER
  48. *          The number of rows of the matrix Q or P**T to be returned.
  49. *          M >= 0.
  50. *
  51. *  N       (input) INTEGER
  52. *          The number of columns of the matrix Q or P**T to be returned.
  53. *          N >= 0.
  54. *          If VECT = 'Q', M >= N >= min(M,K);
  55. *          if VECT = 'P', N >= M >= min(N,K).
  56. *
  57. *  K       (input) INTEGER
  58. *          If VECT = 'Q', the number of columns in the original M-by-K
  59. *          matrix reduced by DGEBRD.
  60. *          If VECT = 'P', the number of rows in the original K-by-N
  61. *          matrix reduced by DGEBRD.
  62. *          K >= 0.
  63. *
  64. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  65. *          On entry, the vectors which define the elementary reflectors,
  66. *          as returned by DGEBRD.
  67. *          On exit, the M-by-N matrix Q or P**T.
  68. *
  69. *  LDA     (input) INTEGER
  70. *          The leading dimension of the array A. LDA >= max(1,M).
  71. *
  72. *  TAU     (input) DOUBLE PRECISION array, dimension
  73. *                                (min(M,K)) if VECT = 'Q'
  74. *                                (min(N,K)) if VECT = 'P'
  75. *          TAU(i) must contain the scalar factor of the elementary
  76. *          reflector H(i) or G(i), which determines Q or P**T, as
  77. *          returned by DGEBRD in its array argument TAUQ or TAUP.
  78. *
  79. *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
  80. *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  81. *
  82. *  LWORK   (input) INTEGER
  83. *          The dimension of the array WORK. LWORK >= max(1,min(M,N)).
  84. *          For optimum performance LWORK >= min(M,N)*NB, where NB
  85. *          is the optimal blocksize.
  86. *
  87. *  INFO    (output) INTEGER
  88. *          = 0:  successful exit
  89. *          < 0:  if INFO = -i, the i-th argument had an illegal value
  90. *
  91. *  =====================================================================
  92. *
  93. *     .. Parameters ..
  94.       DOUBLE PRECISION   ZERO, ONE
  95.       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  96. *     ..
  97. *     .. Local Scalars ..
  98.       LOGICAL            WANTQ
  99.       INTEGER            I, IINFO, J
  100. *     ..
  101. *     .. External Functions ..
  102.       LOGICAL            LSAME
  103.       EXTERNAL           LSAME
  104. *     ..
  105. *     .. External Subroutines ..
  106.       EXTERNAL           DORGLQ, DORGQR, XERBLA
  107. *     ..
  108. *     .. Intrinsic Functions ..
  109.       INTRINSIC          MAX, MIN
  110. *     ..
  111. *     .. Executable Statements ..
  112. *
  113. *     Test the input arguments
  114. *
  115.       INFO = 0
  116.       WANTQ = LSAME( VECT, 'Q' )
  117.       IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'P' ) ) THEN
  118.          INFO = -1
  119.       ELSE IF( M.LT.0 ) THEN
  120.          INFO = -2
  121.       ELSE IF( N.LT.0 .OR. ( WANTQ .AND. ( N.GT.M .OR. N.LT.MIN( M,
  122.      $         K ) ) ) .OR. ( .NOT.WANTQ .AND. ( M.GT.N .OR. M.LT.
  123.      $         MIN( N, K ) ) ) ) THEN
  124.          INFO = -3
  125.       ELSE IF( K.LT.0 ) THEN
  126.          INFO = -4
  127.       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  128.          INFO = -6
  129.       ELSE IF( LWORK.LT.MAX( 1, MIN( M, N ) ) ) THEN
  130.          INFO = -9
  131.       END IF
  132.       IF( INFO.NE.0 ) THEN
  133.          CALL XERBLA( 'DORGBR', -INFO )
  134.          RETURN
  135.       END IF
  136. *
  137. *     Quick return if possible
  138. *
  139.       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
  140.          WORK( 1 ) = 1
  141.          RETURN
  142.       END IF
  143. *
  144.       IF( WANTQ ) THEN
  145. *
  146. *        Form Q, determined by a call to DGEBRD to reduce an m-by-k
  147. *        matrix
  148. *
  149.          IF( M.GE.K ) THEN
  150. *
  151. *           If m >= k, assume m >= n >= k
  152. *
  153.             CALL DORGQR( M, N, K, A, LDA, TAU, WORK, LWORK, IINFO )
  154. *
  155.          ELSE
  156. *
  157. *           If m < k, assume m = n
  158. *
  159. *           Shift the vectors which define the elementary reflectors one
  160. *           column to the right, and set the first row and column of Q
  161. *           to those of the unit matrix
  162. *
  163.             DO 20 J = M, 2, -1
  164.                A( 1, J ) = ZERO
  165.                DO 10 I = J + 1, M
  166.                   A( I, J ) = A( I, J-1 )
  167.    10          CONTINUE
  168.    20       CONTINUE
  169.             A( 1, 1 ) = ONE
  170.             DO 30 I = 2, M
  171.                A( I, 1 ) = ZERO
  172.    30       CONTINUE
  173.             IF( M.GT.1 ) THEN
  174. *
  175. *              Form Q(2:m,2:m)
  176. *
  177.                CALL DORGQR( M-1, M-1, M-1, A( 2, 2 ), LDA, TAU, WORK,
  178.      $                      LWORK, IINFO )
  179.             END IF
  180.          END IF
  181.       ELSE
  182. *
  183. *        Form P', determined by a call to DGEBRD to reduce a k-by-n
  184. *        matrix
  185. *
  186.          IF( K.LT.N ) THEN
  187. *
  188. *           If k < n, assume k <= m <= n
  189. *
  190.             CALL DORGLQ( M, N, K, A, LDA, TAU, WORK, LWORK, IINFO )
  191. *
  192.          ELSE
  193. *
  194. *           If k >= n, assume m = n
  195. *
  196. *           Shift the vectors which define the elementary reflectors one
  197. *           row downward, and set the first row and column of P' to
  198. *           those of the unit matrix
  199. *
  200.             A( 1, 1 ) = ONE
  201.             DO 40 I = 2, N
  202.                A( I, 1 ) = ZERO
  203.    40       CONTINUE
  204.             DO 60 J = 2, N
  205.                DO 50 I = J - 1, 2, -1
  206.                   A( I, J ) = A( I-1, J )
  207.    50          CONTINUE
  208.                A( 1, J ) = ZERO
  209.    60       CONTINUE
  210.             IF( N.GT.1 ) THEN
  211. *
  212. *              Form P'(2:n,2:n)
  213. *
  214.                CALL DORGLQ( N-1, N-1, N-1, A( 2, 2 ), LDA, TAU, WORK,
  215.      $                      LWORK, IINFO )
  216.             END IF
  217.          END IF
  218.       END IF
  219.       RETURN
  220. *
  221. *     End of DORGBR
  222. *
  223.       END
  224.